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DYNAMIC  RESPONSE  OF  A  WETTING  LIQUID  ENCLOSED  IN 
A  ROTATING  TANK  IN  A  ZERO-G  ENVIRONMENT  AND 
SUBJECTED  TO  VARIOUS  DISTURBING  FORCES 

By 

Sui-kwong  Paul  Pao 

April,  1967 

Chairman:   Dr.  Julius  Siekmann 

Major  Department:   Engineering  Science  and  Mechanics 

During  the  past  few  years,  a  growing  interest  in  the  investiga- 
tion of  various  problems  in  the  area  of  the  mechanics  of  contained 
fluids  under  reduced  effective  gravity  has  been  observed,  stimulated 
by  the  practical  needs  of  spacecraft  operations.   Among  others,  the 
question  of  the  dynamic  response  and  stability  of  a  large  body  of  fluid 
in  a  tank  placed  into  subgravity,  is  one  of  considerable  importance. 
The  matter  of  concern  has  to  do  with  the  storage,  transportation,  and 
utilization  of  liquid  propellants  in  the  operation  of  spacecraft  power- 
plants  . 

The  present  study  deals  with  the  exploration  of  the  nature  of 
liquid  interface  oscillations  in  a  rotating  tank  and  the  possibility 
of  rotational  stabilization. 

The  system  under  consideration  consists  of  a  right,  circular, 
cylindrical  tank  of  finite  length,  partly  filled  with  an  Invlscld, 
Incompressible,  wetting  liquid.   The  whole  system  spins  about  the  tank 
axis  with  constant  angular  speed  In  a  near  zero-g  environment. 
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The  large  scale  vapor  bubble  takes  a  stable  equilibrium  shape 
which  is  similar  to  an  elongated  spheroid  and  located  symmetrically 
about  the  axis  of  rotation.  When  perturbing  forces  are  prevailing, 
the  liquid  body  and  the  free  interface  will  oscillate  about  the  equilib- 
rium configuration. 

The  equations  of  motion  are  linearized.  Surface  tension  and 
rotational  velocity  component  are  both  essential  in  the  treatment  of 
the  problem  and  are  taken  fully  into  account. 

The  self-sustained  oscillations  are  governed  by  an  elliptic 
differential  equation  for  the  perturbation  pressure  field.   In  this 
case,  the  bubble  is  approximated  by  a  prolate  spheroid  which  is  embedded 
into  a  spheroidal  coordinate  system.   Hence  associated  Legendre  func- 
tions can  be  employed  for  the  series  expansion  of  the  solution.  A  spe- 
cial method  has  been  employed  in  order  to  account  for  the  homogeneous 
boundary  condition  at  the  walls  of  the  tank.   The  resulting  eigenvalue 
problem  for  the  relative  frequency  of  oscillation  is  nonlinear.   The 
eigen frequencies  are  obtained  through  iteration  procedure.  They  are 
all  real  and  greater  than  two.   The  first  two  eigen frequencies  can  be 
computed  with  great  accuracy.   The  oscillations  are  stable.   The  ampli- 
tude of  oscillation  is  large  near  the  equator  of  the  bubble.   For  very 
slowly  rotating  systems  the  eigenfrequencies  tend  to  accumulate  near 
the  critical  frequency. 

For  a  forced  oscillation,  induced  by  a  reduced  gravity  field 
of  constant  magnitude  and  direction  in  an  inertial  frame  of  reference, 
the  relative  frequency  of  oscillation  is  less  than  two.   The  governing 
differential  equation  is  hyperbolic,  hence  the  method  of  analysis  is 
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different.      Consequently  the   structure   of   the   flow  field   is    different 
from   that  of   the   above   mentioned  elliptic   case.      The  mathematical    prob- 
lem  is    transformed    into    integro-dif ferential   equations   which  can  be 
integrated   numerically  by  means   of   Picard's   method   of   successive   approx- 
imations.     It   is    found   that   a   steady  solution    in  a  cylindrical   tank  of 
finite    length  does   not  exist.      However,    a   solution  exists    in  a   cylin- 
drical  tank  of    infinite    length.      The   perturbations   extend   to   infinity. 
A  numerical   example    for   this    forced   oscillation   problem  has   been   carried 
out. 

We   find    from   this   analysis    that   the  rotation   of   the   system  has 
a   profound    influence   on   the  dynamic   response   of   the    liquid   body.      For 
a   system  without   rotation,    a   small   disturbance    is    liable    to  cause 
instability  of   the   equilibrium  configuration.      The   surface   tension   pro- 
vides  only    limited  means    to  damp   out   any  dynamic   disturbance.      For   a 
system  with   rotation,    all   the  small   oscillations    of   the    interface   are 
stable.      A  small   disturbance  will   be   readily   absorbed   dynamically  with- 
out  destroying   the   stable   configuration.      If   viscous   effects   were    taken 
into  account,    the   oscillations   can   be  effectively  dissipated   through 
the    induced   secondary   flow. 


CHAPTER  I 
INTRODUCTION 


During  the  past  few  years,  a  growing  interest  in  the  investiga- 
tion of  various  problems  in  the  area  of  the  mechanics  of  contained 
fluids  under  reduced  effective  gravity  has  been  observed,  stimulated 
by  the  practical  needs  of  spacecraft  operations.   Among  others,  the 
question  of  the  dynamic  response  and  stability  of  a  large  body  of  fluid 
in  a  tank  placed  into  subgravity,  is  one  of  considerable  importance. 
The  matter  of  concern  has  to  do  with  the  storage,  transportation,  and 
utilization  of  liquid  propellants  in  the  operation  of  spacecraft  power- 
plants  . 

In  low  gravitational  environments,  forces  like  surface  tension 
and  the  centrifugal  force,  induced  by  slow  rotation  of  the  fluid,  will 
have  a  dominating  effect  on  the  large  scale  equilibrium  configuration 
and  the  dynamics  of  a  fluid  system.   For  example,  a  right  circular 
cylindrical  tank,  partially  filled  with  a  wetting  liquid,  spins  about 
its  axis  of  revolution  with  constant  angular  speed  and  is  placed  in 
a  weak  gravitational  (zero-g)  field.   Then  the  surface  tension  at  the 
interface,  together  with  the  centrifugal  force,  causes  the  vapor  cavity 
(bubble)  inside  the  vessel  to  take  an  elongated  spheroid- like  shape, 
situated  symmetrically  about  the  axis  of  rotation.   Equilibrium  con- 
figurations for  various  constant  angular  speeds  have  been  studied  by 


r  -  * 
Rosenthal  11  J.    A  concise  survey  on  the  literature  to  this  field  is 

given  in  a  review  article  by  Habip  [2]. 

The  present  study  deals  with  the  exploration  of  the  nature  of 
liquid  interface  oscillations  in  a  rotating  tank  and  the  possibility 
of  rotational  stabilization. 

The  chosen  frame  of  reference  is  fixed  in  the  tank  and  rotates 
with  it.   The  perturbation  velocities,  the  interface  wave  amplitude, 
and  the  disturbing  forces  are  assumed  to  be  small.   Thus,  in  the  equa- 
tions of  motion,  the  equation  of  continuity,  and  the  boundary  condi- 
tions terms  involving  quadratic  or  higher  orders  of  the  perturbation 
quantities  are  neglected. 

The  fluid  is  assumed  to  be  inviscid  and  incompressible. 
Surface  tension  and  rotational  velocity  components  are  both  essential 
in  the  study.   They  are  taken  fully  into  account  as  far  as  a  linearized 
theory  will  allow. 

For  this  rotating  fluid  oscillation  problem  there  are  two  fre- 
quency ranges  for  which  the  disturbances  are  of  an  entirely  different 
character.   Let  uu  be  the  dimensionless  frequency  of  oscillation  which 
is  measured  relative  to  the  constant  rotation  of  the  tank-liquid  system. 
Then  for  -u  >  2,   the  type  of  the  governing  partial  differential  equa- 
tion is  elliptic,  while  for  uj  <  2,  the  flow  field  is  described  by  a 
hyperbolic  equation.   The  method  of  analysis  and  the  physical  interpre- 
tation of  the  flow  phenomena  are  entirely  different. 


Numbers  in  brackets  refer  to  the  Bibliography. 


For  the  elliptic  case,  the  bubble  is  approximated  by  a  prolate 
spheroid  which  is  embedded  into  a  spheroidal  coordinate  system.  Hence, 
associated  Legendre  functions  can  be  employed  for  the  series  expansion 
of  the  solution.  A  special  method  has  been  employed  in  order  to  account 
for  the  homogeneous  boundary  condition  at  the  wall  of  the  tank.  The 
resulting  eigenvalue  problem  for  the  relative  frequency  of  oscillation 
is  nonlinear,  since  the  eigenfrequency  is  included  aa  a  parameter  in  the 
formulation  of  the  differential  equation  as  well  as  the  boundary  condi- 
tions. The  eigenfrequencies  are  obtained  through  a  successive  approx- 
imation procedure.  This  procedure  proved  to  be  stable  for  the  numerical 
computations  performed  for  a  sequence  of  examples. 

The  result  for  this  case  indicates,  that  for  the  circumferential 
mode  m  =  2,  the  first  eigenfrequency  greater  than  two  corresponds  to 
the  first  mode  of  vibration  in  the  meridian  plane.   Hence,  for  higher 
modes  of  oscillation  the  flow  field  is  certainly  elliptic  in  nature. 

For  a  forced  oscillation,  induced  by  a  reduced  gravity  field  of 
constant  magnitude  and  direction  in  an  inertial  frame  of  reference,  the 
relative  frequency  of  oscillation  is  less  than  two.  The  governing  dif- 
ferential equation  is  hyperbolic,  hence  the  method  of  analysis  is  dif- 
ferent.  Consequently,  the  structure  of  the  flow  field  is  different 
from  that  of  the  above  mentioned  elliptic  case.  The  mathematical  prob- 
lem (Cauchy  and  Goursat  problem)  is  transformed  into  integro-dif ferential 
equations  which  can  be  integrated  numerically  by  means  of  Picard's 
method  of  successive  approximations.   It  is  found  that  a  steady  solution 
in  a  cylindrical  tank  of  finite  length  does  not  exist.   However,  a 


solution  exists  in  a  cylindrical  tank  of  infinite  length.  The  pertur- 
bations extend  to  infinity.  A  numerical  example  for  this  forced  oscil- 
lation problem  has  been  carried  out. 

In  the  present  work  we  consider  only  oscillations  which  are 
symmetric  with  respect  to  the  equatorial  plane  of  the  bubble,  i.e., 
symmetric  with  respect  to  the  variable  £  in  the  cylindrical  coordi- 
nate system,  and  with  respect  to  the  variable  0  in  the  prolate  spher- 
oidal coordinate  system. 


CHAPTER  II 


FORMULATION  Of  THE  PROBLEM 


2 . 1      Derivation   of   the    linearized  equations   of  motion 

The  equations   of  motion  governing   the    fluid   flow   in   a  rotat- 
ing  frame   of   reference   can  be  written    in   the   following   form    [3] 


2t  r       ' 


(1) 


where   q  denotes  the  relative  (perturbation)  velocity  vector,  fl  the 
angular  velocity  vector,  r   the  position  vector  from  the  origin  of  the 
rotating  coordinate  system  to  the  point  occupied  by  the  fluid  particle, 
and   F  the  external  force  per  unit  mass  of  the  fluid.  We  further 
denote  the  density  of  the  liquid  by  o   a"d  the  pressure  by  p.   In  the 
present  problem  the  basic  fluid  flow  is  in  a  state  of  rigid  body  rota- 
tion with  respect  to  the  inertia  frame  of  reference  and  hence  is  at 
rest  with  respect  to  the  rotating,  nontrans lating  frame  as  referred  to 
in  this  investigation.   By  neglecting  the  quadratic  terms  of  the  per- 
turbation velocity  components  a  linearized  system  of  equations  is 
obtained.  This  system  can  be  written  in  a  cylindrical  coordinate 
system  (r,9,z)  as 
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with  fu,0,w]  as  the  components  of  the  perturbation  velocity  vector  and  with 

fF  ,F  ,F  )   as  the  components  of  the  external  force. 
n  9  z 

The  equation  of  continuity  reads 

-L-L(rfl)  +lli  +^  =  o  •  (2) 

The  system  of  differential  equations  can  be  made  dimension- 
less  by  choosing  the  semi-minor  axis  of  the  undisturbed  bubble  radius 
C  as  the  length  scale  and  ,ft    as  the  time  scale.   Further,  the  act- 
ing external  perturbing  forces  in  our  case  are  assumed  to  originate 
from  an  oscillatory  force  field  which  is  transverse  to  the  direction 
of  the  tank  axis.   It  is  assumed  further  that  this  field  has  a  constant 
direction  in  the  inertia  frame  of  reference.   Then  the  pressure  can  be 
written  as 

P  =  ±?n2c2  {t?  +  2v\  +6pc7)e  +Pol>  (3) 

where      6   ,      is    a   constant   and   the   time   dependent,    second   term  of   the 
oL 

right-hand   side    takes   account   of   the   perturbation    force   effects.      The 
equations   of   motion   are   then   reduced   to 


9t  3  7] 

dv  J_3p 
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while  the  equation  of  continuity  becomes 


=  0. 


(5) 


2 .2     Derivation   of   the    field  equation 

By  means   of  simple  elimination,    relations   between  each  veloc- 
ity component   and   the   pressure  may  be   obtained   from  the  equations   of 
motion,   namely, 
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Now  these   relations    can  be   used   to  eliminate   the  velocity  components 
from  the  equation  of  continuity   (5).     Thus,    the  governing  field  equa- 
tion  for  the   fluid   flow  is 


9tz 


(7) 


This  partial  differential  equation  exhibits  some  very  distinct 
properties.   Let  us  assume  that  the  perturbation  pressure  can  be 
written  in  the  form 


y>(r\.9.t,,  t )   =  p  (t),  e  ,t,)e 


(s)t 


(8) 


where,   of   course,    the   real   part  of   the   right-hand  side   has    to  be 
used   for   physical    interpretation. 


With   this   setup   the   field  equation  for   the  perturbation   pres- 
sure   is   transformed   into 

1JL(7}JI)+     I      *l    ,jlz±tl=    o       ,  (9) 

7}  n  l  hi\ '     if  2&2        w2    di2 


and   it   follows   that  the  foregoing  equation   is   of  elliptic   type,  when 
uj  >  2,    and  of  hyperbolic   type,   when     <u  <  2.      The  last   term  on   the    left- 
hand   side   of    (9)   vanishes    if     iju  =   2 .      Similarly,  the   velocity  components 
can  then  be  expressed  explicitly   in  terms  of  the  perturbation  pressure 
and    its   derivatives: 


2  i£     a'<-ut 


(10a) 


(10b) 


(10c) 


Now  let  the  perturbation  velocity  (u,v,w)  be  divided  into  two 
parts  (Uj,  Vj,  Wj)  and  (u2>v2,w2),  such  that 


'?f     'nut 
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The   velocity   component    (u.,v    ,w    )    is    irrotational    because    it    is    t 


he 


gradient   of   a   scalar    function.      For   the    velocity  component    (u   ,v    ,w   ), 
we   have 


curl  V    =    curl   v2   = 
I  \  zALill  2     ft]?**    f     4i      ?2        2       i   tl]t^  [   2   ft]  Art) 

where 

V  =  (  U,V,    W  )   W    v£  =  (  Ui  ,  V2   .  Wij  • 

Hence,  the  velocity  field  can  be  characterized  as  follows; 

a)  For  u)  >  2,  the  magnitude  of  the  rotational  velocity  com- 
ponent is  smaller  than  the  magnitude  of  the  irrotational  one.  Their 
ratio  tends  to  zero  at  the  rate   1/ui  as  u)  tends  to  infinity. 

b)  For  u)  <  2,  the  magnitude  of  the  rotational  velocity  com- 
ponent is  larger  than  the  magnitude  of  the  irrotational  one.   Their 

2 
ratio  tends  to  infinity  of  the  order   0(l/u)  )   as  u>  tends  to  zero. 

The  vorticity  effect  is  so  dominating  that  the  velocity  field  may  have 

to  take  up  a  cellular  structure  in  order  to  make  the  flow  in  the  fluid 

domain  dynamically  possible.      The   hyperbolic   type   of   the   field  equation 

permits   such  discontinuities    in    its   solution.      An  example   can   be   found 

in   Phillips    [U]. 

c)  For  to  =  2,  the  partial  differential  equation  is  indetermi- 
nate in  the  sense  that  the  ^-dependence  is  arbitrary.  Any  solution  of 
the  form 

whereK(T|,9)  satisfies  the  field  equation  and  the  boundary  conditions 
and  any  functions  WAQ)   and  S7  (C)  which  fulfill  the  boundary 
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conditions  are  admissible.  The  solution  is  not  unique.   Hence,  in  the 
linearized  theory  such  a  flow  is  unstable. 

For  the  case  -u  >  2   the  affine  transformation  of  the  space 
coordinates  (Figure  1) 


1    =   V 

e  =  q 


Co 


J  <*l-*t 


1?   > 


(11) 


will    reduce   the   field  equation    (9)    to   the    Laplace   equation 


1  dV  W*V      Vie2      3/i2  ~     ' 


(12) 


\                            \ 

i 

V 

V=Vo 

/ft  =  const.             s' 

0 

\                /  \o(  =  Const . 

/  / 

/     / 

^                      \ 

Figure    1.      A  meridian   section   of    the   deformed   cylindrical 
coordinate   system    (    n   ,    6    ,  ^    )   and    the    prolate 
spheroidal    coordinate    system    (    <<   ,    A    ,  9     ). 
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Another  coordinate  transformation  is  needed  for  the  analysis 
of  this  problem.   The  cylindrical  coordinates  (T],9,  ^,)  are  transformed 
into  the  prolate  spheroidal  coordinates  (or,  P»6)  by  means  of 


1  =f  A^-OO-f)   ' 

/U  =  f  c(jB  ,  <13) 

0  =  9 


where     a     is   a  scale  constant  which  will   be  determined    later. 

In   this    coordinate   system   the   Laplace  equation   for   the   pres- 
sure  field  becomes 

It  should  be  noted  that  the   field  equation   is   separable   in  the  cylin- 
drical  as  well   as    in   the   prolate   spheroidal   coordinate   system. 

2.3     The   boundary  conditions    for   the  elliptic   case 

The  next  step   in  our  analysis    is   to  formulate   the  boundary 

conditions   at   the   rigid   tank  wall   and   the   free    liquid   vapor    interface. 
At   the   plane  end  disks   of  the  tank  the  normal  velocity,   which 

is    in   the    ii-direction,   must  vanish.     Thus 

At  the   cylindrical   portion  of  the  tank  wall  the  normal  velocity  com- 
ponent,  which  is    in  the  Tl-direction,   must  vanish.     Hence, 

97?  7?     29  '         l0  <16> 
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The  bubble  interface  9hall  first  be  specified  in  the  spher- 
oidal coordinate  system  before  the  boundary  conditions  are  set  up. 
The  approximate  bubble  shape  is  assumed  to  be  a  spheroid  which  has  the 
same  ratio  of  the  semimajor  axis  to  the  semiminor  axis  as  the  exact 
shape  of  the  bubble.   The  semiminor  axis  is  normalized  to  one  in  the 
dimensionless  coordinate  systems.   By  taking  into  account  the  stretch- 
ing due  to  the  affine  coordinate  transformation  the  bubble  can  be 
embedded  into  the  spheroidal  coordinate  system  as  one  of  the  coordinate 
surfaces   a  =  o  =  constant,  where   «   is  given  by 

OC0       „     60 


=  K  ,        •  (17) 


The  scale  constant   a  which  appears  in  the  relations  defin- 
ing the  coordinate  transformation  between  (^.G.u)  and  (o,  R,P)  is  then 


a  =  ,        •  (18) 


Hence  the  perturbation  interface  can  be  written  in  the  form 

lUlt 


a. 


m     A.       ♦ 


<*i    (p.  *)«  •  (19) 


2 

where    the    function      o  (P,9)      is   small   such   that      |l  o.     ||  <"<  (a     -    D, 

where       ||  a  denotes    the   norm  of      a,  • 

Now  a   bounding  surface      F(a,8i9,t)    =0     of   a   continuous   medium 
is    a  material   surface.      Mathematically  this   condition   can   be   written   as 

A 


■jf  f  F  (a.jj.e.t)}  -    0    i  (20) 


with 

AE-_lL  +  lIifL    lLAl    l£if.. 
dt       n       2*  dt       -a/3  dt     ?e  dt 


(21) 
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In  the  present  case  we  have 

F(*,j}.e.t)=  oc  -ot0  -  oL,(p.9)eiU)t   =o 


(22) 


hence 


^.  =  ^_lW0c,Me)elWt-^A^t-^i^i-el£dt=o.     (23) 
M         At  '         39  <K       3jB  at 

Neglecting  quadratic  terms  tn  the  above  equation  there  results 

—  -  iudix(&.e)elut  =  o   .  (2U) 

dt        ' 

The  perturbation  velocity  component  in  the  or-direction, 

v   ,  may  be  written  in  terms  of  the  pressure  and  its  derivatives. 
a 

It  turns  out  that 

where  h   denotes  a  scale  factor.   Then  the  above  results  can  be 
Of 

summarized  in  the  following  equation: 

(""O^  -^^^ -77-^90  "i?  '  (26) 

A  dynamical  boundary  condition  shall  also  be  satisfied  besides 
the  above  kinematical  boundary  condition.  The  pressure  at  both  sides 
of  the  interface,  together  with  the  surface  tension  force,  must  be  in 
equilibrium  with  the  centrifugal  and  inertia  forces.   Now,  the  pressure 
in  the  liquid  has  been  defined  previously  by  equation  (3) 
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In  the  cavity,  the  pressure  is  a  constant,   p  „    •  Then  the  boundary 

oo 

condition  at  the  interface  reads 


^iiVf7,2  +  2p)  +  epcT,€^+Jlt) 


I         # 


♦QiL-kWrj  •    (27) 


cl=o<0  +  <x, 


For  the  perturbation  configuration  of  the  bubble  we  have 

I     =-2La(l+o|TT)JK1-l)ri-^     .  <"> 


hence 


^[(^oc.-pK,+^) 


=  fui-0fi-fl)  +£■***,  0-p) 


(29) 


The   following  expression   is   arrived   at  by  substituting 
the    above    quantity    into   the   definition  of    the   pressure    in   the    liquid 
and   by  subtracting   the    quantities   which   are   originally   in  equilibrium 
(static   equilibrium  configuration   of   the   rotating  tank-fluid   system 
with     no  disturbance   forces   acting)   from  both  sides: 


ifflfc'jf'  2<W<-/Jl)  +  2p,j  +  £p^et(ilt  +  el-iTrj'-J>(30) 


We    further  define   the    following  dimens ionless   numbers 
c        P^C*   .  a       -     cfc'         , 

E  =  ~TT  '  B     "     t 

(31) 


where      B      denotes    the   ratio  of    the    relative  magnitude   of    the   centrifugal 
force   to   the    surface    tension   and      B      the    so-called   Bond   number. 
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The   number  E   is   directly  related   to  the   bubble-shaped   factor 

Their  relation   is   shown   in  Figure   2. 

■ 
K 


The  shape  factor  K  versus  the  dimensionless  number  B. 
For  the  limiting  case  E  =  0.5,  the  bubble  becomes 
an  infinite  cylindrical  tube. 


Then,  finally,  the  dynamic  boundary  condition  can  be  written  as 


2 


i(Qt  +  e) 


=  -(J*-J)  ■> 


(32) 


with   pQ  as  the  perturbation  pressure  to  be  evaluated  at  the  interface. 

The  change  of  curvature  J*-  J  is  a  lengthy  expression  in  terms 
of  0,  9,  0^(0,9),  and  the  first  and  second  derivatives  of  a,  (0,9) 
with  respect  to  0.   The  procedure  of  calculation  follows  the  method 
outlined  in  Struik  [5]. 


16 


All  the  energy  transfer  and  other  important  dynamical  effects 
are  contained  in  the  boundary  conditions  at  the  free  interface.  The 
eigenf requencies  are  determined  by  this  boundary  condition. 


CHAPTER   III 

SOLUTION  FOR  THE   ELLIPTIC   CASE 

3.1  Representation  of  the  solution  of  the  field  equation 

Since  the  Laplace  equation  (14)  is  separable  in  the  prolate 
spheroidal  coordinates,  the  perturbation  pressure  may  be  written  as 

!(*,£, *)  =  A(«0M(|9)N(e)-  <33> 

After  a  few  steps  of  standard  computation  the  solution  of  the  field 
equation  is  obtained  in  the  following  forms 

P(*.j».*)»Z£    (  Amn?«(«)  +S™Q~M)tn(f>)el™d  '  (3<o 

In  the  above  representation,  the  functions  ?„<»)»  Pi/fl^' 

and  0m(a)  are  the  Legendre  associated  functions  of  the  first  and  the 
vn 

second  kind,  respectively.   The  coefficients  A^  and  Bmn  are  constants 

to  be  determined  by  the  boundary  conditions. 

im9 
The  dependence  on  the  variable   6  is  given  by  the  terms  e   , 

where  the  number  m  is  an  integer  which  takes  the  values  0,  1,  2 

Physically,  these  are  the  fundamental  modes  as  viewed  on  each  of  the 

circular  cross  sections  transverse  to  the  axis  of  rotation.  The  case 

m  =  0  indicates  that  each  circular  cross  section  of  the  bubble  expands 

or  contracts  radially  with  the  same  amount  for  all   6  .  The  case  m  =  1 
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shows  that  each  circular  cross  section  moves  as  a  whole  transverse  to 
the  axis  of  rotation.   In  this  case  the  dynamic  effects  are  asymmetric. 
If  the  perturbation  pressure  is  integrated  along  any  concentric  circu- 
lar path  inside  the  fluid  domain,  the  result  is  non-zero.   There  is 
a  net  resultant  force  acting  on  the  fluid  body  bounded  by  this  path. 
On  the  other  hand,  if  the  acceleration  field,  induced  by  the  rotation, 
is  symmetric,  the  bubble  is  freely  floating  inside  the  fluid  and  hence 
the  surface  tension  offers  no  resistance  to  an  asymmetric  disturbance, 
except  locally  for  a  small  portion  of  the  entire  interface.   Had  there 
been  no  external  force  field  acting  on  this  system,  there  would  be  no 
restoring  force  to  sustain  such  an  oscillation.   Hence  such  an  oscilla- 
tion is  possible  only  if  some  transverse  inertia  force  field  is  present. 
In  an  inertial  frame  of  reference  a  force  field  realizable  by  this  mech- 
anical system  has  always  a  relative  frequency  of  oscillation  less  than 
two.   Therefore  for  this  mode  of  oscillation  the  flow  field  is  hyperbolic 
in  nature  and  does  not  belong  to  the  present  elliptic  case.   For  the 
case  m  =  2,  the  bubble  undergoes  a  breathing  motion.   For  higher  modes 
the  motion  of  the  fluid  is  of  the  same  nature  as  for  m  =  2. 

3.2  Matching  the  exterior  boundary  condition 

The  solution  of  the  Laplace  equation  as  represented  in  the 
prolate  spheroidal  coordinates  forms  a  complete  system  of  eigenfunc- 
tlons  at  the  bubble  interface   <y  =  a     .   At  the  rigid  cylindrical  tank 
wall,  the  corresponding  sequence  of  functions,  as  given  by  this  repre- 
sentation, is,  however,  not  complete.   It  is  incapable  of  satisfying 
■  boundary  condition  in  general  at  the  exterior  boundary.   Nevertheless, 
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it  is  possible  to  match  a  homogeneous  boundary  condition  at  this  bound- 
ary by  means  of  a  method  described  in  this  section. 

It  may  be  found  also  from  the  computation  of  the  difference  of 
curvature,   J  -  J  ,  that  there  is  no  interaction  between  the  various 
modes  of  e     due  to  the  surface  tension  effect.  Hence,  each  mode  of 
different  values  of  m  can  be  handled  separately.   From  the  discussions 
of  the  previous  section,  only  the  cases  where  m  s  2  are  considered. 
Now  the  summation  with  respect  to  m  may  be  dropped  from  the  represen- 
tation of  the  perturbation  pressure.  Then  the  number  m  remaining  in 
the  expressions  in  question  acts  as  a  parameter  rather  than  an  index. 

In  the  following  a  system  of  eigenf unctions  will  be  obtained  in 
the  cylindrical  coordinate  system  such  that: 

a)  The  eigenfunctions  are  harmonic  inside  the  tank; 

b)  Each  eigenfunction  satisfies  part  of  the  required 
boundary  conditions  at  the  rigid  wall  of  the  tank; 

c)  The  system  of  eigenfunctions  is  complete  over  the 
entire  exterior  boundary,  i.e.,  the  cylindrical  wall 
and  the  plane  end  disks. 

Let  it  be  required  that  the  scalar  product  formed  between  P  and  the 
first  N  eigenfunctions  over  the  entire  surface  of  the  exterior  bound- 
ary are  all  zero,  and  further  let  N  pass  to  infinity.   Since  the  sys- 
tem of  eigenfunctions  is  complete,  the  perturbation  pressure  satisfies 
the  homogeneous  boundary  condition  in  the  mean  at  the  exterior  boundary. 

The  Laplace  equation  as  written  in  the  cylindrical  coordinate 
system  is 

n  ±2L+>*l   n.      .  (35) 
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Let     P     be   represented  as 


I  (-n.s./A)  =   U(t))  V(e)  W(/j) 


(36) 


With  this  scheme  of  separation  of  variables,  the  following  system  of 
uncoupled  ordinary  differential  equations  is  obtained 


ri£d«) -<**?)" 


^dr]  »  '  dT\ 


=    0  •> 


dlV 

de1 

dLW 
dp1 


+    maV    =  0 


+      A2W  =   0. 


(37) 


The  eigenvalues  for  the  parameter  m  are  naturally  determined  as  the 

sequence  of  integers  0,  1,  2,  ...  .   The  eigenvalues  of  the  parameter 

2  2 

X   are  obtained  separately  for  the  case   X  >  0  and  the  case   X   «■"  0. 

These  two  cases  are  different  in  character.  When  the  parameter  m   is 

not  zero,  the  value   X  =  0   is  not  an  eigenvalue  for  this  system  of 

differential  equations  and  boundary  conditions. 


(A)   The  case   X  >  0 

Under  the  requirement  that  the  pressure  field  is  symmetric  with 
respect  to  y,  and  that  the  solution  should  be  regular  inside  the  tank, 
the  solution  of  the  Laplace  equation  (35)  can  be  represented  as 


(38) 


where      I        is    the  modified   Bessel    function  of    the   first    kind   nnd    the 
m 

a    's   are   constants.      The  eigenf unctions   must   fulfill    the    boundary 
n 

condition 
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(  co1-  4)  M   =  i  w  ^   =  °        at   /*  =  * 


i.e.,  at  the  plane  end  disks  of  the  tank,  where  i       --==== 

4  w  -  4 


Hence 


£  a-ni-iUl^Inl  **■<)')  Sin  (*n/u)e  =0       .  (39) 


n=i 


In  order  to  have  a  system  of  nontrivial  solutions,    it  is  necessary 
that 


S*n   (*n/j)  |  .  .,     =   0 


(U0) 


Hence  the  eigenvalues  are  determined  by 

An  =-^7T  n  =  1,2,3,  ...  .     (41) 

This  means   physically  that   the  normal   velocity  at  the  ends  of 
the   tank  vanishes.      The   cylindrical   boundary   is    left  unrestricted   and 
is   free   to  oscillate.     The  normal  velocity  at  the   free  cylindrical 
surface    is   given  by 

(of-A)U  =X    flnO^O  {f  (lm-.(A„T]0)+Inn-ieAnT]o))  + 


This  sequence  of  e i gen funct ions ,  representing  the  normal  veloc- 
ity at  the  entire  exterior  boundary,  has  the  following  properties: 
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a)  The  normal  velocity  at  the  ends  of  the  tank  vanishes; 

b)  The  eigenf unctions  are  orthogonal  to  each  other  over  the 
entire  surface  of  the  exterior  boundary; 

c)  The  sequence  of  eigenfunctions  is  complete  when  its 
domain  is  restricted  to  the  cylindrical  portion  of  the 
exterior  boundary. 

2 
(B)  The  case   X   <  0 

The  main  purpose  for  constructing  these  sequences  of  eigen- 
functions is  that  the  normal  velocity  representation,  derived  from  the 
eigenfunctions  of  the  perturbation  pressure  at  the  entire  exterior 
boundary,  form  a  complete  system.   Hence  the  choice  of  the  boundary 
condition  for  the  eigenfunctions  is  more  or  less  free.   In  the  present 
case  there  is  no  straightforward  boundary  condition  that  can  be  pre- 
scribed such  as  in  case  (A).   The  following  boundary  condition  has  been 
selected  to  start  with 


^7)M  =  Tlo  (U3) 


The  reasons  for  choosing  this  boundary  condition  are  the  following: 

a)  This  is  the  weakest  boundary  condition  which  allows 
a  choice  among  a  variety  of  other  conditions; 

b)  The  eigenfunctions  thus  obtained  do  not  depend  on  the 
parameter  id; 

c)  After  orthogonal izat ion  with  respect  to  the  first 
sequence  of  eigenfunctions  as  obtained  in  the  case  (A), 
each  of  the  eigenfunctions  obtained  in  this  case  will 
satisfy  the  homogeneous  boundary  condition  at.  the 
cylindrical  boundary  in  the  mean; 

d)  The  normal  velocity  representation  at  the  exterior 
boundary,  corresponding  to  the  eigenfunctions ,  forms  a 
sequence  of  functions  which  are  orthogonal  and  complete 
over  the  domain  restricted  to  the  ends  of  the  tank. 
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In  the  present  case  the  representation  of  the  function  P  is 
given  by 

n-i  ' 


where     J       is   the  Bessel  function  of   order     m     and  the     b    's     are 
m  n 

constants.     The  corresponding  velocity  components     u     and     w     are 

Cco2-4)m=X  b„(.'An^7i;^nTj)+~Jffl^7))]cosh^)elM+meJ  (45) 


(W1-4)'M/=X    bnO>n  w)J,n(An^)Sin/)  (  An  u)  e 


(46) 


and 


(^^^[at^^J^A^lsin^Y)^^4^.  (47) 


n  =  l 


In  order  to  satisfy  the  prescribed  boundary  condition  it  is  neces- 
sary to  put 


Hence  the  eigenvalues  are  the  zeroes  of  the  function  J'(*  A  )• 

°  tn  n  o 

Thus  the  two  sequences  of  eigenf unctions  are  completely  defined, 

When  the  proper  values  of  the  coefficients  a   and  b   are  given, 

r   r                             n  n 

these  eigen functions ,  evaluated  at  the  exterior  boundary,  can  be 
written  as: 


2k 


Sequence    (A) 


t  ty  =     COS  (  An  M-} 


*'-l 


Im  UnT)) 


(*n1J.) 


Sequence    (B) 


at  the  cylindrical  wall, 


at  the  plane  end  disks. 


Jm  (Anl.) 


l|>n    =  — — — *-=}    cosh(*nu)    at   the   cylindrical  wall, 


at   the   plane  end   disks. 


WO 


(50) 


with  n  =  1,2,  ...  . 

It  is  possible  to  put  these  two  sequences  of  eigenf unctions 
together  in  order  to  form  a  single  sequence  of  eigenfunctions  such 
that  the  normal  velocity  representation  given  by  these  eigenfunctions 
forms  a  complete  system  of  functions  at  the  exterior  boundary.   Now 
there  arises  the  question  as  to  how  the  functions  in  these  two 
sequences  are  to  be  denumerated  into  a  single  sequence  of  eigenfunc- 
tions.  In  the  present  study  their  relative  importance  is  pointed  out 
by  the  following  consideration. 

By  returning  to  the  prolate  spheroidal  coordinates,  the  zeroes 

of  the  function   Pm(P)   are  uniformly  distributed  along  the  cr-varlable 

n 

where     cp     is   defined   as      cos  a  3  B    .      In   the   three-dimensional   con- 
figuration,   the   perturbation    pressure   would   vanish   along  some    conlr.il 
surface      p   °   const.      Let   the   coordinate   surface    passing   the    intersec- 
tions   of    the   cylinder  with   the   end   disks   be   denoted   by      6   "   +  P 
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For  sufficiently   large  numbers  of     n,    the  number  of  waves  of   the   func- 
tion    P   (8)     outside   the  hyperboloids   and   the  number  of  waves    inside 
n 

approaches   a  constant  ratio     v,   given,  by 


V  - 


^    -  ?°        ,  (51) 


where  cos  cp  =  8  ,   and  co  <  ■» 


3o  "  po'   a"u  ^o  "  2  * 


For  example: 


when  co  =  30°,  then  v  =  2.000. 

o 

Now  the  relative  importance  of  the  boundary  condition  at  the  cylindrical 
wall  and  the  boundary  condition  at  the  end  disks  of  the  bubble  oscilla- 
tion problem  is  assumed  to  follow  the  same  ratio  v.  Then  the  new 
sequence  of  eigenfunctions  is  denumerated  from  the  sequence  (A)  and  the 
sequence  (B)  in  such  a  way  that  in  the  first  N  terms  of  the  new 
sequence,  the  ratio  of  the  members  from  sequence  (A)  to  the  members 
from  the  sequence  (B),  shall  always  be  the  rational  number  closest  to  v. 

For  the  above  mentioned  example,  the  denumeration  of  the  new 
sequence  of  functions   {$  }  is 

$,=  r",   " 

£3=  h        ' 

$4=  ^S   ' 

etc. 

As  said  before,  the  normal  velocity  representation,  given  by 

the  eigenfunctions  of  the  sequence  (A),  forms  a  complete  sequence  of 
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functions  at  the  cylindrical  portion  of  the  boundary  and  vanishes  at 
the  disk  ends,  while  the  normal  velocity  representation,  given  by  the 
eigen functions  of  the  sequence  (B)  forms  a  complete  system  of  functions 
at  the  end  disks  and,  moreover,  is  a  set  of  monotone  functions,  which 
vanish  very  rapidly,  as  n  becomes  large,  over  the  half  lengths  of 
the  cylindrical  portion  of  the  boundary.   In  the  new  sequence   f  $  }, 
the  members  from  the  sequence  (A)  and  the  sequence  (B)  approach  a 
constant  ratio.   It  is  easy  to  verify  that  the  new  sequence  of  eigen- 
f unctions  gives  a  complete  system  of  functions  for  the  normal  veloc- 
ity representation  of  the  exterior  boundary. 

The  eigenf unctions  in  the  sequence  [   $  }  can  be  orthogon- 
al ized  by  using  the  Hilbert-Schmidt  procedure.   In  the  present  study, 
only  the  homogeneous  boundary  condition  has  to  be  satisfied.  The 
sequence   f  $  }  can  be  used  directly.   The  results  obtained  are 
equivalent  to  those  obtained  by  using  the  orthogonal ized  sequence  of 
e  igenfunctions . 

Now  the  equivalent  statement  of  the  homogeneous  boundary  con- 
dition is 

<J>  iif^e^fC-f  e)^-o,  n=  lp  2>  3>  mmm  (52) 

where  E  denotes  the  entire  surface  of  the  tank. 

In  the  above  integrals,  the  function  P  and  the  eigenfunc- 
tions   f  *n  }   are  expressed  analytically  in  two  different  coordinate 
systems.  Hence,  these  integrals  can  be  handled  conveniently  only 
in  a  numerical  fashion. 
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3.3  The  dynamic  (interior)  boundary  condition 

Since  various  modes  of  different  values  of  m  are  not  coupled, 
the  interior  boundary  condition  (32)  for  m  >  2  will  reduce  to 

gE  {$*.*,(  l-f*)*?.)   =  -Cj"-J)  •  (53) 

It  should  be  noted  that  the  term   B  ">]  €         in  equation  (32) 
belongs  to  the  case  m  ■  1.  In  equation  (53)  the  perturbation  pres- 
sure at  the  interface  p   ,  the  difference  of  curvatures,   J  -  J,  and 

o 

the  perturbation  of  the  bubble  shape   ©(..(0,9),  can  all  be  expressed  in 
terms  of  a  series  of  the  associated  Legendre  functions  with  unknown 
constant  coefficients.  The  constant  uj  is  the  eigenfrequency  para- 
meter.  For   p   we  find 

?.  -  £  (*„„£«  ♦  B-.<j;  M)Ol»  e''""4^  .        cm 

The  dynamic  boundary  condition  is  evaluated  at  the  original  equilibrium 
boundary  of  the  bubble,  a   =  a  ,    thus  this  boundary  condition  is,  apart 
from  the  factor  el<,«t"fmfl>|  a  function  of   0  only.  After  some  lengthy 
computations,  for  the  quantities   J*-  J,  oC^fl)  ,     '  ;   ,    and  -   U-  t 
there  results 

*(<-/>'Xi-^)f-.-')-7^-J    + 

t(«'.""'!7'»,,y)f^' !  • 
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«.^=7fe)Lf*^  (56) 


with 


a'0l)  =  -Ok  -0  ^o2-o2    . 

fly  =  -«<WH)r^i)Vyflc-i)^.iX5<H)-2^iXi*i> 
a»>  =  - 2o*0v<-.)  ♦fAw-iWwf-0-2f/w-i)rw-() » 


(?) 


2n   _,* 
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a'?   =  2(4-  "?<  +  0  tJ^W  +  O-fl  +  ftJW+Z)-"2   • 

4- 

a<5)  =  n2  +   7"    +  12        * 

a r';  =  "  2  f  n  -  m  +  / )  ol]  Uo2  +  2)    • 

f      _  Ki^  -I 

7C  o(. 

^    =  4  1     K2t«.2-0   "  11 


The  factor  (1  -  B2)  vanishes  at  (9  *  *  1.  However,  the  terms 

involving  this  factor  in  the  denominator  remain  regular  due  to  the  fact 

that  the  associated  functions  Pm(B)  have  a  zero  of  the  same  or  higher 

n 

order  at  these  points  when  m  s  2.   For  the  same  reason  it  is  permis- 
sible to  multiply  the  equation  of  the  boundary  condition  by  a  factor 
(a2  -  (32)2'(1  -  P2).  By  doing  so,  the  boundary  condition  is  reduced 
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to  a  summation  of  terms  of  P  (8)  with  polynomials  of  B  as  the  coef- 
ficients.  It  appears  in  the  form 

where  J>     r  (8),   i  =  1,2, 3,  4,   are  polynomials  of  8.   These  poly- 
nomials are  obtained  by  substituting  equations  (54),  (55),  and  (56) 

into  equation  (53)  and  multiplying  the  entire  expression  by  a  factor 

2    2  2       2 
(aQ   -  S  )  •(!  -  p  )•  The  highest  order  terms  in  the  polynomials  are 

8 
B  .  ThLS  boundary  condition  can  be  reduced  further  to  a  series  of 

Pn(0)  with  constant  coefficients  by  means  of  the  recurrence  formula 

P'"^  =^br{^-w+')C^  *r-m)/£p)j .      (60) 

Then,    finally,    the  dynamic  boundary  condition  can  be  written    in  the 
form 


oo      K=ntJ 
n  =  m    K=n-8  "  * 


<P  =  0     »  (61) 


where    the      f        (or   ;m,nj),    and      e        (n   ;m,w),    are   constants   determined 
by      n>  ,   m,   n,    and   containing     <u      linearly  as   an  eigenparnmeter.      The 

coefficients   of    the   polynomials    "t>  (ft),    as   well   as    the   constants 

J     n 

f„  (or  ;m,j;)  and  g  (a  ;m,'u),  which  follow  from  equation  (59)  by 
repeated  application  of  equation  (60),  are  obtained  numerically  for 
given   values   of     m,    a  ,      and 


They   are   generated   by   the   computer   program    and  will   not   be 
reproduced  here. 
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A  system  of  linear  algebraic  equations  in  terms  of  A    and 

mn 

B    can  be  obtained  by  forming  the  scalar  products  of  this  boundary 
mn 

condition  and  the  eigenfunctions   P  (B)   for  n  £  m.   Together  with 
the  exterior  boundary  condition,  there  are  sufficient  algebraic  equa- 
tions to  determine  the  eigenvalues  of  u>,  and  subsequently  the  eigen- 

solutions  of  A   and  B   ,  which  determine  the  mode  shape  of  the  oscil- 
mn      mn 

lations  of  the  bubble  interface. 

In  concluding  this  section,  there  is  one  remark  which  should 
be  made  concerning  the  effects  of  the  rotation  of  the  entire  system 
and  the  surface  tension  of  the  oscillations  of  the  bubble.   In  the  final 
form  of  the  interior  boundary  condition  there  is  a  second  summation  sign 
which  sums  over  k  for  nine  terms.    In  an  elementary  oscillation  prob- 
lem, the  second  summation  sign  is  not  present.   Each  simple  mode  of  the 
harmonic  oscillations,  as  given  by  the  eigenfunctions,  is  distinct, 
such  as  in  the  case  of  the  vibrations  of  a  string  with  fixed  end  points. 
In  the  present  linearized  theory  of  oscillations,  the  surface  tension 
effect  induces  an  interaction  among  nine  simple  modes  of  harmonic  oscil- 
lations, while  the  rotational  effect  alone  will  induce  an  interaction 
among  three  consecutive  simple  modes  only. 

3.4  Numerical  solution 

If  the  series  representation  of  the  perturbation  pressure  is 
truncated  to  n   terms  only,  the  boundary  conditions  can  be  reduced  to 
a  system  of  linear  algebraic  equations  in  terms  of  the  unknown  constant 


* 
In  the  present  case  the  oscillations  are  symmetric  with  respect 

to  0.  Hence  only  the  eigenfunctions  with  even  values  of  (n-m)  have  non- 
zero coefficients.   Only  these  terms  are  counted. 
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coefficients,  and  the  eigenparameter  m  is  contained  in  these  equa- 
tions. Since  the  oscillations  are  symmetric,  only  the  terms  with  even 
number  of  (n-m)  are  non-zero.  There  are  only  (n-m)  ♦  2  unknown  con- 
stants, A   and  B   .   Each  of  the  boundary  conditions  will  render 
mn      mn 

—(n-m)  ♦  1   linear  equations  of  the  unknown  constants.  The  system  of 
linear  equations  is  determinate. 

The  reduction  of  the  boundary  condition  to  a  system  of  alge- 
braic equations  is  given  in  the  following. 

The  exterior  boundary  condition  (52)  yields 

Z    f  9*k    Arvk       +4*k    Bmk]     =  0  >      L    =     1,2,       ,j(n-m)  +  l     (62) 
K  =  m 


where 

while  the  interior  boundary  condition  (61)  becomes 

Z  {0*oi'n.w)Amk  *1l°(«..m.u)hnk]*°*     i-*.«*2,-.«*«  (63) 

where   f„   (cy  ;m,fu)   and   g_  (a   ;m,u))  are  the  same  as  before. 
K    o  TQ    o 

Furthermore,  the  coefficients  can  be  written  as 

(.)  #  en   j_  , •  tn      . 

j-  (<*0.m.  u)   ■  a,k   -  w  olk    , 


g'j;,(otclm.w)  =  at.. 
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At  this  point  it  is  convenient  to  introduce  matrix  notation  to 
handle  the  system  of  algebraic  equations.   All  the  matrices,  G,  H,  A    , 
A(2),  B(1),  B(2),  defined  below,  are   fy(n-m)  ♦  l]  X  t^0-™)  *  l1 
square  matrices.  The  vectors  X  and  Y  are  two  column  matrices  with 
— (n-m)  ♦  1  elements.   They  are 


(,\         l      CO)  CO  #(|) 

A0)=ja..j5    a,y  =  a*       , 


CO  ,*Cl) 


ttco       Juc,)l  W      -  b 


X    =  (  *i)  ,      xt  =  A 


m<<  ■» 


A».{a?J.  a"-^.        y    -  f  ttl '      *  =  B"*  ' 


(64) 


1  ,  2,   .   .   .  ,   {(m-n-)+  1  , 

m    +■   2(  i  -0         ? 


=    m 


+  2  Cj  -  0 


Then  the  exterior  boundary  condition  together  with  the  interior  boundary 
condition  can  be  written  in  the  following  concise  form: 


I 


A'°  !    A(z> 


H 


l_ 

u 


gCrtj   gCZ) 


0 


I 


0 


X 

y 


=  0 


(65) 


This  matrix  equation  can  be  solved,  employing  standard  nota- 
tion as  follows : 

First  we  have 

GX  ♦  HY  =  0,  (66) 


34 


which  yields 


Further 


X  =   -  G~   HY.  (67) 


(A(1)    -IB(1))X-    (A(2)-Ib(2))Y=0,  (68) 

en  (u 


upon  substituting  for  X,  we  have 


or 


(A(1)  -iB(1))(-G-1HY)  ♦  (A(2)  -IB(2))Y-0, 


f(A(2)-A(1)G-1H)  -  -   (B(2)-BG-1H)-|Y  =  0.  (69) 


Then  with  the   definitions 


there   results 


R  =  A^-A^V'h   and   S  »  B(2).B(1)o\ 


(R  -  -  S)Y  =  0.  (70) 


Equation  (70)  can  be  written  as 


(R  -  -  S)S"lSY  =  0, 


thus,    finally 


(RS_I    -  -  E)SY  =    0,  (71) 

U) 

where   E    is    the   unit  matrix. 

The   eigenvalues     —     and   the  eigenvectors      Z   =    (SY)      can  be 
found    from  the    last  equation.      Then   the  eigenvectors   can  be   expressed 
in   terras   of    the   original      X  and   Y: 

Z   =   SY, 
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yields 


and  hence  by  equation   (67) 


Y  =  S_1Z,  (72) 


X  =  -  g"1hy. 


For  any  specific  example,  where  dimensions  of  the  tank,  the 
density  of  the  liquid,  the  constant  speed  of  rotation,  and  hence  the 
static  equilibrium  shape  of  the  bubble  are  given,  numerical  results 
can  easily  be  obtained  by  performing  the  computations  according  to  this 
theory  on  a  digital  computer. 

In  the  numerical  computation  the  scale  of  the  affine  trans- 
formation 


J  a>*-4  Ln    r  j~^~a 

has  to  be  chosen  at  the  beginning.   Then  according  to  the  obtained 
eigenvalues  of  uu  this  scale  of  the  affine  transformation  can  be  com- 
puted employing  the  above  given  formula  for  p,.  The  new  value  of  this 
9cale  factor  is  then  used  for  the  computation  of  a  more  accurate  value 
of  id. 

It  is  important  to  point  out  that  a  different  scale  factor 
exists  for  each  different  eigenvalue  of  to*   For  very  large  values  of  to 
the  scale  factor  is  close  to  1. 
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3 .5  Results  and  discussion 

Numerical  computat ions  have  been  made  for  the  case  m  =  2, 
which  is  the  lowest  circumferential  mode  under  the  scope  of  the  ellip- 
tic case.   The  result  shows  that  the  first  eigenfrequency  above  2.000 
corresponds  always  to  the  first  mode  of  oscillation  in  the  meridian 
plane.   For  any  specific  example  we  are  able  to  compute  the  first  and 
the  second  eigenf requencies  accurately.  We  may  conclude  that  all  the 
natural  frequencies  of  the  bubble  oscillations  are  greater  than  two 
and  their  corresponding  perturbation  pressure  field  is  elliptic  in 
nature,  since  we  have  obtained  the  lower  bound  of  these  frequencies. 

Some  of  the  numerical  results  are  presented  in  Table  1,  Fig- 
ure 3,  and  Figure  4.   The  computation  is  performed  for  a  sequence  of 
bubble  shapes  while  the  tank  dimensions  remain  the  same.   The  first  and 
the  second  eigenf requencies  for  each  shape  factor  K  are  given  numer- 
ically in  Table  1.   They  are  also  plotted  versus  the  shape  factor  K  in 
Figure  3.   Typical  configurations  for  each  mode  in  a  meridian  section 
are  given  in  Figure  4.   It  is  interesting  to  note  that  the  excitation 
of  the  interface  wave  is  much  more  pronounced  in  the  neighborhood  of 
the  equator  than  in  the  pole  regions  of  the  bubble. 

In  the  numerical  computations   n   is  taken  as  large  as  m  ♦  12 
such  that  up  to  seven  non-zero  terms  in  the  series  are  included.  At 
least  five  non-zero  terms  shall  be  taken  in  order  to  obtain  an  accurate 
result  for  the  eigenf requenc ies .   This  is  not  surprising  because  the 
interactions  among  the  fundamental  harmonic  functions  are  very  strong. 

The  third  eigenfrequency  may  be  obtained  also.  However,  the 
accuracy  is  considerably  less  than  the  first  two,  due  to  the  limited 
capacity  of  the  computer. 
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TABLE  1 


THE  FIRST  AND  SECOND  EIGENFREQUENCIES  FOR  BUBBLES 
WITH  DIFFERENT  SHAPE  FACTORS  IN  A  CYLINDRICAL 
TANK  WITH  FIXED  DIMENSIONS 


The  shape  factor       The  first         The  second 

K  eigenfrequency     eigenfrequency 


1.100 
1.200 
1.400 
1.600 
1.800 
2.000 
2.200 


2.296 

2.409 

2.401 

2.558 

2.516 

2.780 

2.620 

2.977 

2.754 

3.235 

2.941 

3.605 

3.168 

4.007 

,  v     =   2.000, 


o 


The   dimensions   of   the   tank:      -j 

L      t   =   3.000. 
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l_  /.i                                                                         -i 

3.50 

It) 

71      =   2.000 

'o 

I      =   3.000 

m     =   2 

3.00 

2.50 

^         ^ 

s. 

^ First  eigenfrequency 

/ 

/ 

-    Second  eigenfrequency 

/ 

2.00 

1.00  1.20 


2.00 


K 


Figure    3.      The    first   and   second   e igenf requeue len    plott.    I 
versus    the    bubble    shnpt     (actor    K. 
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3.0 

2.0 

1.0 

P 

0 

0.5 

1.0 

1.0 

' 

Figure  4.   The  typical  mode  shapes  in  the  meridian  plane 

corresponding  to  the  first  and  the  second  eigen- 
frequencies  of  the  oscillations  of  the  bubble 
about  its  position  of  stable  equilibrium. 


CHAPTER    IV 
SOLUTION  TO  THE   HYPERBOLIC  CASE 

4. 1      Introduction 

In   this   chapter   the   dynamical   response   of   the   rotating 
fluid   system  with  respect   to  the   perturbation   due   to  a  constant 
reduced   gravity   field   transverse   to  the   tank  axis  will   be   studied. 
From  the   discussion   on   page    18,   we   know  that   for   the  mode   of  oscil- 
lation    m   =    1,    the  motion  has   to  be   accompanied   by  an  external 
force   field.      The   governing  equation    is   hyperbolic.      Hence   the 
method  of   analysis    and   the  mathematical    formulation   of   the   problem 
are   different   from  the   elliptic   case   studied   above. 

For   the   physical   conditions   given   above,    the    frequency  of 
the    perturbation      o  relative   to  the   rotating  system   is   one.      How- 
ever,   the   analysis    for   this    particular   case,    as   persued    in   the    fol- 
lowing sections   of   this   chapter,    is   valid   also   for   any   frequency 
in   the   range      0  <r  ■»  <   2 . 

The   mathematical    formulation  here    is   written    in   cylin- 
drical  coordinates   and   the    "characteristic   coordinates."     The    latter 
will   be   defined    in   this   chapter. 

Now  the   side   conditions    for   a  well-posed   problem   in   the 
sense   of   Hadamnrd   are   different    for  hyperbolic   and   elliptic  equa- 
tions.     For   the   present   problem,   we   have   to   prescribe    kinematic 
boundary  conditions   nt   the    cylindrical    wall,    the   disk  ends   and   the 
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interface,  as  well  as  a  dynamic  condition  at  the  interface.   Such 
a  system  of  boundary  conditions  is  too  stringent  for  the  solution  to 
exist.  We  shall  see  in  section  4.4  that  the  boundary  condition  at 
the  ends  of  the  tank  has  to  be  relaxed.   Hence  we  are  considering 
a  cylinder  of  infinite  length  instead  of  one  with  finite  length. 
Since  the  differential  equation  governing  the  pressure 
field  is  hyperbolic,  there  are  real  characteristic  surfaces  in  the 
flow  domain.   Across  these  surfaces  the  normal  derivatives  may  be 
discontinuous.   Hence  the  velocity  field  will  suffer  a  finite  jump 
at  the  same  location.   The  work  by  Oser  [6]  provides  a  good  example. 

4.2  Formulation  of  the  problem 

The  field  equation  in  cylindrical  coordinates  (T],0,C)  for 
the  case   m  =  1   and   uu  =  1   is  given  by 

fl  +  ±JL  -- Lf-3^  =0    ■ 

?f     1  37)      t)2         at,1 

The  first  order  derivative  term  in  this  equation  can  be  eliminated 
by  means  of  the  transformation 

Fdj.eX)  =  7~'  iCn-O.O-  (73) 

Then   the   partial   differential  equation   for      $      takes   the   form 

tk     _3Jl     =    A±9      .  (74) 

37)2        3  3CZ         4   f  * 

This  transformation  has  been  carried  out  for  the  sake  of  conven- 
ience regarding  later  numerical  computations. 
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For   the   above   hyperbolic   equation    (74),    there   are   two 
families   of   real   characteristic   curves    through  each   point   of   the 
C^tO   plane.      Both  of   these   families   are   straight    lines. 

These   characteristics   may  be   chosen   as   our   coordinate 
system.      The   characteristic   transformation   in   question   is 

gjj   +(    =  2  fit 

h  -rj  -   i,    =  2/3  T       »  (75) 

where   x  =  const,  and   t  =  const,   are  the  characteristic  lines 
from  each  of  the  two  families. 

Along  the  characteristic  directions,  the  transformation 
does  not  determine  the  second  derivatives  uniquely.   Hence  across 
a  characteristic  line,  the  normal  derivative  may  suffer  a  jump. 
However,  the  function   $   is  supposed  to  be  continuous  for  a  hydro- 
dynamical  problem. 

Through  the  transformation  (75),  the  hyperbolic  differ- 
ential equation  (74)  turns  into  the  normal  form: 

For  this  equation  two  types  of  problems  can  be  posed. 

(A)  The  Cauchy  problem 

The  initial  conditions  are  prescribed  on  a  curve  PQ  which 
is  nowhere  tangent  to  the  characteristic  directions.   Thus  the  values 
of  the  function   *   and  its  normal  derivative  on  this  curve  are 
given.   A  solution  exists  and  is  uniquely  determined  in  the  triangu- 
lar domain   PQR  (Pigure  5). 
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Figure   5.      A  sketch  for   the  Cauchy  problem. 

The  entire   system  of  differential  equations   and   the    initial   conditions 
can  be   combined    into  the   following   integro-dif ferential   equation 
(Garabedian    [7]): 


i(8)  = 7—  +T 


h 


3$ 


ra« 


t+u+t)2 


d%dx 


(77) 


This  equation  can  be  solved  by  an  iteration  technique.   This  process 
is  proved  to  be  convergent.  The  solution  exists  in  the  large  when  the 
equation  is  linear,  as  is  the  present  case. 
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(B)  The  Goursat   problem 

We  may  also   prescribe   one   boundary  condition   on  the   ordinary 
(non-characteristic)   curve   and   one   on   a  characteristic    line.      If   the 
value   of      §      is    prescribed   on   these   curves,    the   problem   is   called   the 
Goursat   problem.      A   solution  exists   also   and    is   uniquely  determined    in 
the   triangular   domain   OTQ   (Figure   6). 


X  =  %(z) 


JZj  (*o,X0) 


Figure   6.      A  sketch   for   the   Goursat   problem, 


This    problem  can  be   formulated   as   an   integral   equation,    namely, 


l(k)  ■  $(?)*l(Q)-l(S)    + 


3$ 


KQK    f(t**)' 


djidt 


(78) 
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This   equation   can  be   solved   by  means   of   Ptcard's   method  of   successive 
approximations   as    in    (A).     When  the   boundary  condition   on   the   ordinary 
curve   is   replaced  by  a  mixed  type  condition 

H    +acx)-^    +  b(x)$  =   0  ,  (79) 

the  existence   and  uniqueness   of   the   solution   remain  unaltered.      The 
problem   is  well    posed.      However,    the   integro-differential   equation   is 
somewhat  complicated.      Let  us   denote   the  ordinary  curve   by 

The   function     \(.t)      is   strictly  monotonic.     Then  the    integro-differ- 
ential equation   in   question   can  be  written   in   the   form: 

ilk)  -  i(Q)  -  |(5)  -  e"Ift#1  [  bft)fW  +a(i$]  •  eIMdx  ,  &     -rg^todz  i 


(80) 
I(x)  =       b(x)4t 

The  solution  of  the  Goursat  problem  for  a  linear  equation  exists  in 
the  large,  as  in  case  (A). 

Numerical  computations  have  been  performed  according  to  the 
above  given  integro-differential  equations. 

4.3  The  boundary  conditions 

The  boundary  conditions  prescribed  for  this  problem  have  been 
derived  in  Chapter  II.   They  are  given  by  equations  (15),  (16),  (24), 
and  (32).  They  are  listed  below  aa  rewritten  in  the  cylindrical 
coordinates  (T\,d,Q). 
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(A)  The  kinematic  boundary  conditions  at  the  walls  of  the  tank: 

!i  +  .Li  *  =  0  at  the  cylindrical  wall, 

(81) 

~ 77"  —  0  at  the  plane  disk  ends. 

(B)  The  kinematic  boundary  condition  at  the  interface: 

Vn     -  S     =  O        '  (82) 


The  normal  velocity  v   and  the  normal  displacement   6  of  the  liquid- 
vapor  interface  are  actually  90  out  of  phase  when   0  and   t  are 
taken  into  account. 


(C)  The  dynamic  boundary  condition  at  the  interface: 

8E(.T\Sn   +p0)  +Bt/  =  -Cj*-J)  • 

where      6_      is   the   component   of      6      in   the   Tl-direction.      This   boundary 
T] 

condition,    especially   its   right  hand   side,    is   very  complicated.     We 
shall   replace    it  by  two  asymptotic   representations.      This   approximation 
produces   some    quantitative   errors,   while    it   demonstrates,    however,    the 
qualitative   natures   of   the   problem  much  more   clearly. 

For   a   small   displacement      6      of   the    bubble    interface,    the   change 
of   curvature   may  be   approximated   by 

r>     t  dlf> 

J    "J=     d?       ' 

where  ds   is  a  length  element  measured  along  the  generator  of  the 
bubble  interface. 
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In  the  neighborhood  of  the  equator  of  the  bubble,      d£  ~  ds . 
Hence  we   have   the   following  asymptotic   representation  of   the   boundary 
condition: 


SET  (  Vn     +  Po)    +    Bt)  =  -— i 


(83) 


In  the  neighborhood  of  the  poles  of  the  bubble  the  curvature 
of  the  bubble  is  pronounced.   The  geometry  in  this  region  is  best 
described  through  the  spheroidal  coordinates  defined  in  the  previous 
chapters  (Figure  7). 


0-60- 


Figure    7.      The   coordinate   system  employed  at   the 
pole   region  of   the   bubble. 


U8 


The  formulas  applicable  to  the  present  situation  are  summarized  below: 


-  "/[« 


fi     -  cos  (f  » 

7)   =  sin  <f  ■> 

t  =  Kcoscf 

d(f  1  A 

-[J=  /J  I  +  (K2-0  5in>        ' 

Near  the  pole  the  value  of  cp   is  small.   The  curvature  may  be  written  as: 

d:6         d%  (  df\z 


The     T|-component   of   the   displacement      6,    6_,      is   also  very  small    in 
this   region.      It   can  be   neglected.      Hence  we  have   another  asymptotic 
representation  for  the  dynamic   boundary  condition: 


d2S    I  d<f  V2 
8£j>c    +   &S/H9    =   _   d(f2V  dsj     *  (84) 


The   quantity      (",    )      approaches   one   as     cp     approaches   zero. 

QS 

The   combination   of   the    kinematic   condition    (82)   and    the   dynam- 
ical  condition    (83)   at    the    interface   provides   us   only  with  one   boundary 
condition   for   the    field  equation    (76). 

(D)  The   condition   of   symmetry  at   the  equatorial    plane 
We   assume   that   the    flow  field    is   symmetric  with   respect   to   the 
equatorial    plane      Q  =   0.      Hence  we   have    the   condition 


11 


■  0 
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There  are  two  more  conditions  for  the  determination  of  the  flow 
field.  These  conditions  will  be  introduced  in  the  next  section  in 
connection  with  the  construction  of  the  solution  to  this  problem. 

4.4  Construction  of  the  solution 

For  a  hyperbolic  equation  discontinuities  of  the  normal  deriv- 
atives are  admissible  across  a  characteristic  line.  We  shall  construct 
a  solution  to  this  problem  with  all  the  possible  discontinuities  in 
mind.   Figure  8  shows  the  division  of  the  flow  field  into  regions 
separated  by  characteristics  where  discontinuities  may  occur.   In  this 
figure  the  point  G  is  located  such  that  the  bubble  is  tangent  to  the 
characteristic  line. 

As  indicated  in  4.3,  we  have  to  postulate  two  additional  condi- 
tions to  construct  the  flow  field. 

a)  In  the  immediate  neighborhood  of  the  equatorial  plane,  the 
flow  field  is  symmetric.   The  perturbation  velocity  normal  to  this 
plane  vanishes.   The  bubble  is  tangent  to  a  cylinder.   Hence  the  per- 
turbation pressure  distribution  given  in  Phillips  [4]  is  exactly  the 
same  as  the  perturbation  pressure  distribution  at  the  equatorial  plane. 
In  our  notation  there  follows 

*l<;  =  o   "  "1  8E  ^  I  If   X}  (85) 

b)  The  perturbation  velocity  at  the  axis  of  rotation  is  finite. 

From  the  expression  of  the  perturbation  velocity  (6),  we  arrive  at 

the  condition 

£  =  0   • 
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Pigure   8.      Division   of   the    flow   field    Into  regions    separated   by 
characteristic    lines  where   discontinuities  may  occur. 
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We  can  now  proceed  with  the  construction  of  the  solution. 

The  condition  4.3(D)  together  with  condition  4.4(a)  form  a  set 
of  initial  conditions  for  the  determination  of  the  flow  in  the  region  I. 
Consequently,  the  generalized  Goursat  problems  can  be  defined  in  the 
regions  II  through  V  successively.  Their  solutions  are  obtained  by 
means  of  equation  (78)  or  equation  (80).   In  the  region  VI,  the  solu- 
tion is  obtained  as  follows. 

The  displacement  of  the  bubble  6  at  cp  =  0  vanishes  for  the 
geometrical  compatibility  requirement.   In  a  neighborhood  of  the  pole 
such  that  T)  <  |  —  |  ,  both  the  inertia  and  the  centrifugal  field 
forces  are  smaller  than  the  surface  tension  and  the  reduced  forces  by 
at  least  two  orders  of  magnitude.  Hence  from  equation  (84)  we  have 


or, 


Al    _  _  R  I"  |  +  (  k2  -  I )  5  in(f  1  sin  if 
d<fz       l 


hence 


i    _  R  f  (  I  +!(K2-0)  s,nc?  +  T^f]   +  Bo-o^  ' 

0  L      J  (86) 


where  the  constant  e       is  determined  by  the  condition  that  the  dis- 

o 

placement  at  G  is  continuous.   By  means  of  equations  (82)  and  (84), 
the  solution  in  region  VI  is  determined. 

When  we  proceed  to  construct  the  solution  in  further  regions, 
we  find  that  the  problem  is  overdetermined  when  the  boundary  condition 
at  the  disk  end  of  the  tank  is  prescribed.  Consequently  a  solution  to 
the  problem  would  not  exist  in  the  entire  flow  domain  at  all.  This 
boundary  has  to  be  removed  in  order  that  a  solution  may  exist .  Once  this 
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boundary  is  removed,  we  see  immediately  that  the  solution  extends  to 
infinity.  Physically  it  means  that  a  steady  state  solution  does  not 
exist  for  a  tank  of  finite  length.  The  solution  exists  for  a  cylin- 
drical tank  of  infinite  length.  In  this  case  the  perturbation  prop- 
agates to  infinity  immediately.  Such  an  unexpected  phenomenon  was  found 
in  a  similar  problem  investigated  by  Benjamin  and  Barnard  [8]. 

The  solution  in  the  regions  VII,  VIII,  IX,  ...  ,  can  easily 
be  determined  by  a  sequence  of  generalized  Goursat  problems.  The  per- 
turbation field  approaches  zero  as  we  extend  the  solution  to  infinity 
along  the   ^-direction. 

A  numerical  example  is  given  in  Figures  9  and  10. 

U. 5  Discuss  ion 

From  the  results  obtained  from  this  chapter,  we  have  the  follow- 
ing physical  picture  for  the  hyperbolic  case. 

The  liquid-vapor  system  is  stable  with  respect  to  the  forced 
perturbation,  since  its  dynamic  response  to  the  perturbation  is  an 
oscillation  about  the  stationary  equilibrium  configuration  with  a 
small  amplitude. 

In  the  region  RSGA  (Figure  10),  there  is  a  strong  exchange 
between  the  kinetic  energy,  the  energy  of  the  perturbation  pressure 
field,  and  the  centrifugal  field  effects.   The  effect  of  surface  ten- 
sion iB  negligible  in  this  region.   There  is  a  strong  secondary  circu- 
lation induced  by  the  disturbing  forces. 

The  disturbances  in  the  region  below  AG  are  small.   Centrif- 
ugal and  inertia  effects  are  negligible  in  this  region.   The  surface 
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tension  provides  the  necessary  adjustment  to  absorb  the  pressure  field 
perturbations.   In  this  region,  the  liquid  particle  performs  mainly 
an  oscillatory  motion  in  the  neighborhood  of  its  own  equilibrium 
position. 


CHAPTER  V 
CONCLUSIONS 

In  the  previous  chapters  we  have  determined  the  dynamic 
response  of  the  rotating  fluid  system  for  the  entire  range  of  fre- 
quencies.  In  conclusion,  we  shall  discuss  here  some  of  the  phys- 
ical implications  of  the  results  obtained  above. 

The  alow  rotation  of  the  system  with  a  constant  angular  speed 
has  a  profound  effect  on  the  dynamic  response.  For  a  rotating  system, 
all  the  small  oscillations  of  the  liquid-vapor  interface,  or  rather, 
of  the  entire  liquid  body,  are  stable.  A  small  transverse  disturb- 
ance to  the  system  will  induce  one  or  several  modes  of  oscillation 
about  the  stable  equilibrium  configuration. 

On  the  other  hand,  under  the  influence  of  a  transverse  constant 
force  field,  the  stability  of  the  configuration  of  the  system  in  the 
neighborhood  of  the  equator  of  the  bubble  is  ensured  by  the  centrif- 
ugal pressure  field  and  the  inertia  force  produced  by  a  small  pertur- 
bation to  the  constant  rotating  base  flow.  The  surface  tension  effect 
is  negligible  in  this  region.   For  a  system  without  rotation,  a  dis- 
turbance containing  such  a  force  component  is  liable  to  excite 
instability. 

Furthermore,  for  a  real  system  viscosity  effects  are  always 
present.   For  a  system  with  rotation,  a  perturbation  will  induce  some 
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secondary  circulation.  The  disturbance  will  be  dissipated  by  means  of 
the  viscous  mechanism.   The  surface  tension  effect  alone  is  a  two- 
dimensional  mechanism.   It  can  restore  a  disturbed  system  to  an  equi- 
librium configuration  in  a  much  longer  period  of  time. 

Finally,  we  shall  present  here  an  example  to  show  the  actual 
magnitude  of  various  quantities  characterizing  a  rotating  tank-liquid 
system  in  low  gravity  environments. 


Example 


Given: 

The  radius  of  the  tank  R  =  100  cm   , 

o 

Semiminor   axis   of   the 

undisturbed  bubble  C     =     40  cm     . 

Bubble   shape   factor  K     =2.00 

The   dimens ion  less   number  E 

corresponding  to  the  bubble 

shape   factor  K  E   =   0.43  •> 

3 
Density  of   the    liquid  p    =   0.070  g/cm 

Surface   tension  T   =   2.00  dyne/cm 


The    liquid   chosen    in   this   example    is    liquid  hydrogen.      From   the   above 
given   data,   we   may  compute   the   corresponding  constant   angular   speed   of 
rotation     n     and   the   allowable  magnitude   of   the    transverse   reduced 
gravity   field: 


SET  /      faH  "20  =  J_    ra<jl  , 


D     <<    SZ 


Assume  that 


Then 


B     € 


6  =  O  \o    x  12  x  C   =  y^-    /^c 


For  this  small  acceleration  the  system  can  cover  a  distance  of  100 
feet  in  fifteen  minutes.  This  would  allow  the  system  to  perform  a 
slow  maneuver  like  changing  the  orientation  of  the  tank,  etc. 
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APPENDIXES 


APPENDIX  A 


EXPRESSION  OF  THE   PERTURBATION   VELOCITY  COMPONENT  v      IN  THE 

a 

PROLATE  SPHEROIDAL  COORDINATES  (a, 0,8) 


The  expression  of  v   may  be  derived  from  the  known  expressions 
of  u  and  w   in  the  deformed  cylindrical  coordinates  (7),  9,^), 


d?  1   1? 


(Al) 


by  means   of   the   defined  coordinate   transformation  relation,    equation    (13). 
Let  us    define 


9* 

9j3 

Vfo.fi 

97, 

97, 

3ty/i) 

3<x 

3/1 

V 

9/ 

According  to  the   coordinate   transformation   defined    in  equation    (13), 
we    find 


9(*.f) 
9^./*) 


~f 


PJl-i82       Jrf'-I 


(A2) 


The  laws  of  transformation  for  the  contravariant  vector 


.  dii       dii  .  .  .  ,  2?       dP  s 

(  — L  )  — *=-  )   and   the   covanant  vector    ( ■  )   are   given   as 

flit       dt  ?T}      d" 


61 


62 


it 
dt 


and 


97) 

2£ 

(A3) 


(AH) 


where   the   superscript     T      in  equation   (A3)   denotes    the  transposed 
matrix,    and 


*n 


m  '-  M   '   Tt 


■      =    TAT     ,         <*>U  -T-      =     VaL       • 


2-  J  ot2--  I     dt 


(A5) 


Hence,    by  using  equations    (Al),    (A2),    (A3),    and    (AH), 
we   obtain 


id         |_    I"     Ku    (dl-\)    3f        _S* 3f_ 

"dt   =     w2-4  L     a1     (dx-fL)  9*       aV-f)  IB 


(A6) 


The   expression   for  V^    follows    immediately   from  equations    (A5)   and    (A6), 


APPENDIX  B 

FLOW  CHART  OF  THE  COMPUTING  PROGRAM  EMPLOYED 
IN  CHAPTER  III 


(  START  } 


READ:   Tank  dimensions. 

Bubble  shape  factor. 

Assumed  scale  of  transformation 


"4 


w*-4  * 


Dimension  statement. 

Equivalence  statement. 

Define  complex  variable  names, 


I 


FUNCTION: 

The  relation  of 

the  shape  factor 

K  versus  the 

dimension  less 

number  E. 

Arithmetic  statement  for  COSH(X). 


Compute  cK0    according  to  formula  (17). 
Determine  constants  of  the  coordinate 
transformation  from  C7),©.^)  to  (ot.ft,8  ), 


I 


Compute  all  constants  related  to  equations  (55), 
(56),  (57),  and  (58).  Store  results  by  means  of 
subscript  variables.  


I 


Construct  a  polynomial  multiplication  segment. 
Compute  P'^(B)         according  to  equation  (55), 


I 


Construct  a  segment  of  program  for  the  repeated 
application  of  the  recurrence  formula  (60). 

Compute    j^  (oL.  ;  m,  U))        and  q  £  (<x0  ,  m  ,  a>) 


I 


Store  the  obtained  / 'Kn)  (<*„ ; m ,  u)  and  $e*(ct.,m.  u)  in 
proper  locations  of  the  matrices  A(10  x  20)  and 
B(10  x  20). 


<D 


<D 
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© 


FUNCTIONS: 
The   associated 
Legendre   functions. 


READ:      The   size 
of   increment   for 
the  numerical 
integration   along 
the  walls   of   the 
tank. 

"l 


FUNCTION:  The 
hypergeometr ic 
function. 


COMPUTE:      The   values   of 

f>>)P;y>and    Q:(«)P?<p 

at  points  on  the  tank  walls. 
About  25  to  50  points  are 
taken  for  a  quarter  of  the 
tank  along  a  meridian  section 
1  ine . 


FUNCTION:   The 
Gamma  function, 


FUNCTIONS:   The  Bessel 
function  J  and  the 
modified  Bessel 
function  I. 


Construct  the  sequence  of  eigen- 
functions  $; .   Compute  their  values 
at  the  assigned  points  on  the  tank 
walls.   Store  result  in  matrix  form. 


Computation  of  the  numerical 
integrations  according  to 
the  formulas 


READ:   Zeroes  of 
the  Bessel  func- 
tion J'.   The 
denumeration  of 
the  sequence  of 
eigen functions . 


Store  the  obtained  values 
of  <j"k  and  -ft\K  in  proper 
locations  of  a  matrix 
ATK(10  x  20). 


© 
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READ:   The  range 
of  variation  of  n. 
Say,  4  to  7. 


Transcribe  n  x  n  blocks 
from  the  matrices  A,  B,  and 
ATK  to  form  the  new  matrices 

(1)    (2)    (1)    (2) 

A    ,  A    ,  B    ,  B    ,  G, 

and  H   in   equation    (65). 


READ:      A   scale 
factor   for  over- 
flow control  . 


Solve   the  matrix  equation    (65) 
by   the   steps    (66)   through   (71), 


SUBROUTINE: 
Matrix   inversion. 


WRITE:   The 

e igenf re quenc tes 


SUBROUTINE:   For  finding 
the  eigenvalues  and  eigen- 
vectors for  a  real  matrix. 


The  eigenvectors  given 
in  terms  of  the  trans- 
formed bas  is . 


\  WRITE : 
\  eigenve 

The    / 

ctors.  J 

1 

(^  Eb 

«0 

Transform  the  eigenvectors 
back  to  the  original  basis 
by  equations  (67)  and  (72), 
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